###### Models and Results ##### #---------------------------------------------#

# This script runs all regression models in the paper and the Appendix, as well as
# figures based on these models.

# load final data set
load("data/data_final.Rda")


##### Table 1: Summary table #####
data_final_fe %>% 
  mutate(tenure_lag_log = log(tenure_lag),
         pop_lag_log = log(pop_lag),
         gdp_pc_lag_log = log(gdp_pc_lag)) %>% 
  select(any_event, all_events,
         v2x_civlib_lag_inv,
         e_polity2_lag,
         elections,
         tenure_lag_log,
         v2x_ex_party_lag,
         v2x_ex_military_lag,
         pop_lag_log,
         gdp_pc_lag_log,
         Capacity_lag,
         diplo_rep_lag) %>% 
  st(out = "latex", file = "tables/table1.tex")


##### Table 2: Regression results #####

# Model 1, domestic repression and fixed effects only
m1_fe_civ <- glm(any_event ~  I(v2x_civlib_lag_inv*10) + 
                   iso3c + year, data = data_final_fe,
                 binomial(link = "logit"),
                 intercept = F,
                 method = "brglmFit", type = "MPL_Jeffreys")
#summary(m1_fe_civ)

# Model 2, adding control variables
m2_fe_civ <- glm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                   e_polity2_lag + elections + I(log(tenure_lag)) +
                   v2x_ex_military_lag  + 
                   v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) +
                   iso3c + year, data = data_final_fe,
                 binomial(link = "logit"),
                 intercept = F,
                 method = "brglmFit", type = "AS_mixed")

# Model 3,interaction with state capacity
m3_fe_civ <- glm(any_event ~  I(v2x_civlib_lag_inv*10)*I(Capacity_lag*10) +
                   e_polity2_lag + elections + I(log(tenure_lag)) +
                   v2x_ex_military_lag  + 
                   v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) +
                   iso3c + year, data = data_final_fe,
                 binomial(link = "logit"),
                 intercept = F,
                 method = "brglmFit", type = "AS_mixed")


# Model 4, adding interaction with diplomatic representation
m4_fe_civ <- glm(any_event ~  I(v2x_civlib_lag_inv*10)*I(diplo_rep_lag/10) +
                   e_polity2_lag + elections + I(log(tenure_lag)) +
                   v2x_ex_military_lag  + 
                   v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) +
                   iso3c + year, data = data_final_fe,
                 binomial(link = "logit"),
                 intercept = F,
                 method = "brglmFit", type = "AS_mixed")

texreg(file = "tables/table2.tex",
       list(m1_fe_civ, m2_fe_civ, m3_fe_civ, m4_fe_civ),
       stars = c(0.001, 0.01, 0.05,0.1),
       booktabs = TRUE,
       label = "main_reg",
       omit.coef = "(year|iso3c|Intercept)",
       float.pos = "h",
       caption = "Main regression results. Relationship between civil liberties (V-Dem) and transnational repression events",
       custom.coef.names = c("Domestic repression",
                             "Polity score", "Elections",
                             "Leader tenure (log)", 
                             "Military dimension index",
                             "Party dimension index",
                             "Population (log)",
                             "GDP per capita (log)",
                             "State capacity index",
                             "Repression x State capacity",
                             "Diplomatic representation abroad",
                             "Repression x Diplomatic representation"))


##### Table A1: fixest package ######

# Model 1, fixed effects only
m1_fixest <- feglm(any_event ~ I(v2x_civlib_lag_inv*10) | 
                     iso3c + year, data = data_final_fe,
                   binomial(link = "logit"),
                   panel.id = ~ iso3c + year,
                   vcov = "NW")


m2_fixest <- feglm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                     e_polity2_lag + elections + I(log(tenure_lag)) +
                     v2x_ex_military_lag  + 
                     v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) |
                     iso3c + year, data = data_final_fe,
                   binomial(link = "logit"),
                   panel.id = ~ iso3c + year,
                   vcov = "NW")

m3_fixest <- feglm(any_event ~  I(v2x_civlib_lag_inv*10)*I(Capacity_lag*10) +
                     e_polity2_lag + elections + I(log(tenure_lag)) +
                     v2x_ex_military_lag  + 
                     v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag)  |
                     iso3c + year, data = data_final_fe,
                   binomial(link = "logit"),
                   panel.id = ~ iso3c + year,
                   vcov = "NW")

m4_fixest <- feglm(any_event ~  I(v2x_civlib_lag_inv*10)*I(diplo_rep_lag/10) +
                     e_polity2_lag + elections + I(log(tenure_lag)) +
                     v2x_ex_military_lag  + 
                     v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10)  |
                     iso3c + year, data = data_final_fe,
                   binomial(link = "logit"),
                   panel.id = ~ iso3c + year,
                   vcov = "NW")


texreg(
  list(m1_fixest, m2_fixest, m3_fixest, m4_fixest),
  file = "tables/table_A1.tex",
  stars = c(0.001, 0.01, 0.05,0.1),
  booktabs = TRUE,
  label = "main_reg",
  float.pos = "h",
  caption = "Robustness tests. Relationship between dome and transnational repression events",
  custom.coef.names = c("Domestic repression",
                        "Polity score", "Elections",
                        "Leader tenure (log)", 
                        "Military dimension index",
                        "Party dimension index",
                        "Population (log)",
                        "GDP per capita (log)",
                        "State capacity index",
                        "Repression x State capacity",
                        "Diplomatic representation",
                        "Repression x Diplomatic representation"))

##### Table A2: Alternative model specifications ######

# Model 2 feglm
m2_fixest_glm <- feglm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                         e_polity2_lag + elections + I(log(tenure_lag)) +
                         v2x_ex_military_lag  + 
                         v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) |
                         iso3c + year, data = data_final_fe,
                       binomial(link = "logit"),
                       panel.id = c("iso3c", "year"),
                       vcov = "NW")

m2_fixest_pois <- fepois(all_events ~ I(v2x_civlib_lag_inv*10)  +
                           e_polity2_lag + elections + I(log(tenure_lag)) +
                           v2x_ex_military_lag  + 
                           v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) |
                           iso3c + year, data = data_final_fe,
                         panel.id = c("iso3c", "year"),
                         vcov = "NW")

m2_fixest_negbin <- fenegbin(all_events ~ I(v2x_civlib_lag_inv*10)  +
                               e_polity2_lag + elections + I(log(tenure_lag)) +
                               v2x_ex_military_lag  + 
                               v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) |
                               iso3c + year, data = data_final_fe,
                             panel.id = c("iso3c", "year"),
                             vcov = "NW")

m2_fixest_ols <- feols(all_events ~ I(v2x_civlib_lag_inv*10)  +
                         e_polity2_lag + elections + I(log(tenure_lag)) +
                         v2x_ex_military_lag  + 
                         v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) |
                         iso3c + year, data = data_final_fe,
                       panel.id = c("iso3c", "year"),
                       vcov = "NW")

texreg(file = "tables/table_A2.tex",
       list(m2_fixest_glm,m2_fixest_pois, m2_fixest_negbin, m2_fixest_ols),
       stars = c(0.001, 0.01, 0.05,0.1),
       booktabs = TRUE,
       label = "main_reg",
       omit.coef = "(year|iso3c)",
       float.pos = "h",
       caption = "Robustness tests. Main Model (Model 2) with different model specifications.",
       custom.coef.names = c("Domestic repression",
                             "Polity score", "Elections",
                             "Leader tenure (log)", 
                             "Military dimension index",
                             "Party dimension index",
                             "Population (log)",
                             "GDP per capita (log)",
                             "State capacity index",
                             "Theta"))

##### Table A3: Stepwise inclusion of controls #####

# Model 2, adding control variables
m1_contr <- glm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                  e_polity2_lag +
                  iso3c + year, data = data_final_fe,
                binomial(link = "logit"),
                intercept = F,
                method = "brglmFit", type = "AS_mixed")

m2_contr <- glm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                  e_polity2_lag + elections +
                  iso3c + year, data = data_final_fe,
                binomial(link = "logit"),
                intercept = F,
                method = "brglmFit", type = "AS_mixed")

m3_contr <- glm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                  e_polity2_lag + elections +  I(log(tenure_lag)) +
                  iso3c + year, data = data_final_fe,
                binomial(link = "logit"),
                intercept = F,
                method = "brglmFit", type = "AS_mixed")

m4_contr <- glm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                  e_polity2_lag + elections +  I(log(tenure_lag)) +
                  v2x_ex_military_lag  + 
                  v2x_ex_party_lag  + 
                  iso3c + year, data = data_final_fe,
                binomial(link = "logit"),
                intercept = F,
                method = "brglmFit", type = "AS_mixed")

m5_contr <- glm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                  e_polity2_lag + elections +  I(log(tenure_lag)) +
                  v2x_ex_military_lag  + 
                  v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) +
                  iso3c + year, data = data_final_fe,
                binomial(link = "logit"),
                intercept = F,
                method = "brglmFit", type = "AS_mixed")

m6_contr <- glm(any_event ~ I(v2x_civlib_lag_inv*10)  +
                  e_polity2_lag + elections +  I(log(tenure_lag)) +
                  v2x_ex_military_lag  + 
                  v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) +
                  I(Capacity_lag*10) +
                  iso3c + year, data = data_final_fe,
                binomial(link = "logit"),
                intercept = F,
                method = "brglmFit", type = "AS_mixed")

texreg(file = "tables/table_A3.tex",
       list(m1_contr, m2_contr, m3_contr, m4_contr, m5_contr, m6_contr),
       stars = c(0.001, 0.01, 0.05,0.1),
       booktabs = TRUE,
       label = "main_reg",
       omit.coef = "(year|iso3c|Intercept)",
       float.pos = "h",
       caption = "Main regression results. Relationship between civil liberties (V-Dem) and transnational repression events",
       custom.coef.names = c("Domestic repression",
                             "Polity score", "Elections",
                             "Leader tenure (log)", 
                             "Military dimension index",
                             "Party dimension index",
                             "Population (log)",
                             "GDP per capita (log)",
                             "State capacity index"))


##### Figure 4: Bivariate relationship ####

ggplot(data_final_fe, aes(x = civlib_inv_diff2, y = any_event)) +
#  geom_rug( sides="b") +
  geom_point(aes(color = factor(any_event)), alpha = .35) +
  xlab("Change in domestic repression (t-1)") +
  ylab("Probability of transnational repression (t)") +
  scale_color_grey() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"),
              formula= y ~ x, se = T,lwd=0.5, colour="black") +
  theme_bw() +
  theme(legend.position = "none")

ggsave("figures/Figure_4.pdf", width = 6, height = 5, dpi = "retina")





#### Figure 5: Predicted probabilities ####

# Based on Model 2

# Calculating predicted values, varying domestic repression and keeping controls at representative values
pre1 <-predictions(m2_fe_civ, by = "v2x_civlib_lag_inv",
                   newdata = datagrid(v2x_civlib_lag_inv = seq(min(data_final_fe$v2x_civlib_lag_inv),
                                                               max(data_final_fe$v2x_civlib_lag_inv),.05),
                                      FUN_factor = unique, FUN_numeric = median),
                   type="response")

# Plotting predicted values
ggplot(pre1) +
  geom_ribbon(aes(v2x_civlib_lag_inv, ymin = conf.low, ymax = conf.high), alpha = .2) +
  geom_line(aes(v2x_civlib_lag_inv, estimate)) +
  scale_x_continuous("Domestic repression (t-1)", expand = c(0,0)) +
  scale_y_continuous("Predicted probability of TR (t)", expand = c(0,0)) +
  coord_cartesian(ylim = c(0, 1), xlim = c(0, 1)) +
  theme_bw()

ggsave("figures/Figure_5.pdf", width = 7, height = 5, dpi = "retina")

#### Figure 6: Predicted probabilities and interactions ####

# Based on Model 3

# Calculating predicted values, varying domestic repression and keeping controls at representative values
pre2 <-predictions(m3_fe_civ, by = c("v2x_civlib_lag_inv", "Capacity_lag"),
                   variables = list("Capacity_lag" = "threenum"),
                   newdata = datagrid(v2x_civlib_lag_inv = seq(min(data_final_fe$v2x_civlib_lag_inv),
                                                               max(data_final_fe$v2x_civlib_lag_inv),.05),
                                      FUN_factor = unique, FUN_numeric = median),
                   type="response") %>% 
  mutate(Capacity_lag = round(Capacity_lag, 2))


# Plotting predicted values
p2 <- ggplot(pre2 %>% filter(Capacity_lag !=-0.06), aes(x = v2x_civlib_lag_inv,
                                                  y = estimate, 
                                                  fill = factor(Capacity_lag))) +
  geom_line(aes(linetype = factor(Capacity_lag))) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2) +
  scale_fill_manual("State capcacity",
                    values = c("grey75", "grey10"),
                    labels = c("mean - 1 sd", "mean + 1 sd")) +
  scale_linetype_manual("State capcacity",
                        values = c(2, 1),
                        labels = c("mean - 1 sd", "mean + 1 sd")) +
  scale_x_continuous("Domestic repression (t-1)", expand = c(0,0)) +
  scale_y_continuous("Predicted probability of TR (t)", expand = c(0,0)) +
  coord_cartesian(ylim = c(0, 1), xlim = c(0, 1)) +
  theme_bw() +
  theme(legend.position = "bottom")

# Based on Model 4

# Calculating predicted values, varying domestic repression and keeping controls at representative values
pre3 <-predictions(m4_fe_civ, by = c("v2x_civlib_lag_inv", "diplo_rep_lag"),
                   variables = list("diplo_rep_lag" = "threenum"),
                   newdata = datagrid(v2x_civlib_lag_inv = seq(min(data_final_fe$v2x_civlib_lag_inv),
                                                               max(data_final_fe$v2x_civlib_lag_inv),.05),
                                      FUN_factor = unique, FUN_numeric = median),
                   type="response")  %>% 
  mutate(diplo_rep_lag = round(diplo_rep_lag, 2))

# Plotting predicted values
p3 <- ggplot(pre3 %>% filter(diplo_rep_lag != 48.96), aes(x = v2x_civlib_lag_inv,
                                                    y = estimate,
                 fill = factor(diplo_rep_lag))) +
  geom_line(aes(linetype = factor(diplo_rep_lag))) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2) +
  scale_fill_manual("Diplomatic representation",
                    values = c("grey75", "grey10"),
                    labels = c("mean - 1 sd", "mean + 1 sd")) +
  scale_linetype_manual("Diplomatic representation",
                    values = c(2, 1),
                    labels = c("mean - 1 sd", "mean + 1 sd")) +
  scale_x_continuous("Domestic repression (t-1)", expand = c(0,0)) +
  scale_y_continuous("Predicted probability of TR (t)", expand = c(0,0)) +
  theme_bw() +
  coord_cartesian(ylim = c(0, 1), xlim = c(0, 1)) +
  theme(legend.position = "bottom")

ggpubr::ggarrange(p2,p3)

ggsave("figures/Figure_6.pdf", width = 9.5, height = 6, dpi = "retina")

##### Figure A3: Alternative dependent variables #####

# Run models for different targets
m2_fe_civ <- glm(any_event ~ I(v2x_civlib_lag_inv*10) + e_polity2_lag + elections +
                   I(log(tenure_lag)) +
                   v2x_ex_military_lag  + 
                   v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) +
                   iso3c + year - 1, data = data_final_fe,
                 binomial(link = "logit"),
                 method = "brglmFit", type = "MPL_Jeffreys")

m_fe_dv_alt1 <- update(m2_fe_civ, I(citizen > 0) ~ .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1) %>% 
  mutate(outcome = "Citizens")

m_fe_dv_alt2 <- update(m2_fe_civ, I(activist > 0) ~ .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1) %>% 
  mutate(outcome = "Activists")

m_fe_dv_alt3 <- update(m2_fe_civ, I(former_government_official > 0) ~ .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1) %>% 
  mutate(outcome = "Former gov. official")

m_fe_dv_alt4 <- update(m2_fe_civ, I(journalist > 0) ~ .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1) %>% 
  mutate(outcome = "Journalists")

m_fe_dv_alt5 <- update(m2_fe_civ, I(opposition > 0) ~ .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1) %>% 
  mutate(outcome = "Opposition")

m_fe_dv_alt6 <- update(m2_fe_civ, I(executed > 0) ~ .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1) %>% 
  mutate(outcome = "Carried out")

m_fe_dv_alt7 <- update(m2_fe_civ, I(threats > 0) ~ .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1) %>% 
  mutate(outcome = "Threatened")

m_fe_dv_alt8 <- update(m2_fe_civ, I(attempts > 0) ~ .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1) %>% 
  mutate(outcome = "Attempted")

# Combine results
results_dv <- mget(ls(pattern="^m_fe_dv")) %>%
  bind_rows() %>% 
  mutate(type = case_when(outcome %in% c("Activists", "Citizens", "Former gov. official",
                                         "Journalists", "Opposition") ~ "Target",
                          outcome %in% c("Attempted", "Carried out",
                                         "Threatened") ~ "Action",
                          T~ NA_character_))

# PLot estimastes
ggplot(results_dv, aes(x = reorder(outcome, -estimate), group = type, y = estimate,
                       color = type, shape = type)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_point(size = 2) + 
  geom_linerange(aes(
    ymin = conf.low,
    ymax = conf.high),
    lwd = .5) +
  #  ggtitle("1. Politicization effects") +
  scale_x_discrete("Estimate") +
  scale_color_manual("",values = c("#E69F00", "#009E73")) +
  scale_shape_discrete("") +
  scale_y_continuous("", limits = c(-.5,3)) +
  coord_flip() +
  theme_bw() +
  theme(legend.position = "right") 

ggsave("figures/Figure_A3.pdf",
       width = 7, height = 5, dpi = "retina")

##### Figure A4: Component indicators #####

# Rescale component indicators
data_final_fe <- data_final_fe %>% 
  mutate(v2csreprss = I(v2csreprss * -1),
         v2csrlgrep = I(v2csrlgrep * -1),
         v2mecenefm = I(v2mecenefm * -1),
         v2meharjrn = I(v2meharjrn * -1),
         v2meslfcen = I(v2meslfcen * -1),
         v2psbars = I(v2psbars * -1),
         v2psparban = I(v2psparban * -1)
)

m0_fe_alt <- glm(any_event ~     +
                   e_polity2_lag + elections + I(log(tenure_lag)) +
                   v2x_ex_military_lag  + 
                   v2x_ex_party_lag  +  log(pop_lag) + log(gdp_pc_lag) + I(Capacity_lag*10) +
                                iso3c + year - 1, data = data_final_fe,
                              binomial(link = "logit"),
                              method = "brglmFit", type = "MPL_Jeffreys")

m_fe_alt1 <- update(m0_fe_alt, . ~ + v2clkill + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt2 <- update(m0_fe_alt, . ~ + v2cltort + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt3 <- update(m0_fe_alt, . ~ + v2mecenefm + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt4 <- update(m0_fe_alt, . ~ + v2meharjrn + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt5 <- update(m0_fe_alt, . ~ + v2meslfcen + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt6 <- update(m0_fe_alt, . ~ + v2cldiscm + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt7 <- update(m0_fe_alt, . ~ + v2cldiscw + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt8 <- update(m0_fe_alt, . ~ + v2psparban + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt9 <- update(m0_fe_alt, . ~ + v2psbars + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt10 <- update(m0_fe_alt, . ~ + v2psoppaut + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt11 <- update(m0_fe_alt, . ~ + v2cseeorgs + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt12 <- update(m0_fe_alt, . ~ + v2csreprss + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt13 <- update(m0_fe_alt, . ~ + v2clslavem + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt14 <- update(m0_fe_alt, . ~ + v2clslavef + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt15 <- update(m0_fe_alt, . ~ + v2clprptym + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt16 <- update(m0_fe_alt, . ~ + v2clprptyw + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt17 <- update(m0_fe_alt, . ~ + v2clfmove + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt18 <- update(m0_fe_alt, . ~ + v2cldmovem + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt19 <- update(m0_fe_alt, . ~ + v2cldmovew + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt20 <- update(m0_fe_alt, . ~ + v2clrelig + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

m_fe_alt21 <- update(m0_fe_alt, . ~ + v2csrlgrep + .) %>% 
  tidy(conf.int = T, conf.level = 0.95) %>% 
  head(1)

# Combine results
results_iv <- mget(ls(pattern="^m_fe")) %>%
  bind_rows() %>% 
  filter(term != "I(v2x_civlib_lag_inv * 10)")

# Add variables names from codebook
vdem_codeb <- vdemdata::codebook %>% 
  select(term = clean_tag, name)

results_names <- left_join(results_iv, vdem_codeb) %>% 
  mutate(index = case_when(term %in% c("v2clkill", "v2cltort") ~ "Physical violence index",
                           term %in% c("v2mecenefm", "v2meharjrn", "v2meslfcen",
                                       "v2cldiscm", "v2cldiscw", "v2psparban", "v2psbars",
                                       "v2psoppaut", "v2cseeorgs", "v2csreprss") ~
                             "Political civil liberties index",
                           term %in% c("v2clslavem", "v2clslavef", "v2clprptym", "v2clprptyw",
                                       "v2clfmove", "v2cldmovem", "v2cldmovew", "v2clrelig",
                                       "v2csrlgrep") ~ "Private civil liberties index",
                           T ~ NA_character_))

#Plot results
ggplot(results_names, aes(x = reorder(name, estimate), y = estimate, 
                          color = index, shape = index)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_point(size = 2) + 
  geom_linerange(aes(
                     ymin = conf.low,
                     ymax = conf.high),
                 lwd = .5) +
  scale_x_discrete("Indicator") +
  scale_color_manual("",values = c("#E69F00", "#009E73",
                                "#0072B2")) +
  scale_shape_discrete("") +
  scale_y_continuous("") +
  coord_flip() +
  theme_bw() +
  theme(legend.position = "bottom")  +
  guides(colour = guide_legend(nrow = 2))
 # facet_grid(~ outcome ) 

ggsave("figures/Figure_A4.pdf",
       width = 7, height = 5, dpi = "retina")


##### Figure A5: Fixed effects #####
fixedEffects = fixef(m2_fixest)
summary(fixedEffects)

pdf("figures/Figure_A5.pdf", width = 9, height = 5)
plot(fixedEffects)
dev.off()


